{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# install.packages('DirichletReg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "heading_collapsed": "false"
   },
   "outputs": [],
   "source": [
    "library(Seurat)\n",
    "library(RColorBrewer) #for brewer.pal\n",
    "library(Matrix) #for Matrix\n",
    "library(DirichletReg)\n",
    "library(data.table)\n",
    "library(tidyverse)\n",
    "library(cowplot)\n",
    "\n",
    "## this function is extracted from analysis.r \n",
    "dirichlet_regression = function(counts, covariates, formula){  \n",
    "  # Dirichlet multinomial regression to detect changes in cell frequencies\n",
    "  # formula is not quoted, example: counts ~ condition\n",
    "  # counts is a [samples x cell types] matrix\n",
    "  # covariates holds additional data to use in the regression\n",
    "  #\n",
    "  # Example:\n",
    "  # counts = do.call(cbind, tapply(seur@data.info$orig.ident, seur@ident, table))\n",
    "  # covariates = data.frame(condition=gsub('[12].*', '', rownames(counts)))\n",
    "  # res = dirichlet_regression(counts, covariates, counts ~ condition)\n",
    "  \n",
    "  #ep.pvals = dirichlet_regression(counts=ep.freq, covariates=ep.cov, formula=counts ~ condition)$pvals\n",
    "\n",
    "  # Calculate regression\n",
    "  counts = as.data.frame(counts)\n",
    "  counts$counts = DR_data(counts)\n",
    "  data = cbind(counts, covariates)\n",
    "  fit = DirichReg(counts ~ condition, data) \n",
    "  \n",
    "  # Get p-values\n",
    "  u = summary(fit)\n",
    "  #compared with healthy condition, 15 vars. noninflame and inflame, 30pvalues\n",
    "  pvals = u$coef.mat[grep('Intercept', rownames(u$coef.mat), invert=T), 4] \n",
    "  v = names(pvals)\n",
    "  pvals = matrix(pvals, ncol=length(u$varnames))\n",
    "  rownames(pvals) = gsub('condition', '', v[1:nrow(pvals)])\n",
    "  colnames(pvals) = u$varnames\n",
    "  fit$pvals = pvals\n",
    "  \n",
    "  fit\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "freq=read.csv('celltype.csv',row.names=1)\n",
    "freq1=as.matrix(as.data.frame(lapply(freq, as.double),row.names=row.names(freq)))\n",
    "cov=read.csv('conv.csv',row.names=1)\n",
    "cov1 = data.frame(condition=factor(cov[rownames(freq),1], levels=c('healthy control','mild(moderate)','severe','convalescence','influenza')), \n",
    "                  row.names=rownames(freq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in DR_data(counts):\n",
      "“not all rows sum up to 1 => normalization forced\n",
      "  some entries are 0 or 1 => transformation forced”\n"
     ]
    }
   ],
   "source": [
    "pvals = dirichlet_regression(counts=freq1, covariates=cov1, formula=counts ~ condition)$pvals\n",
    "colnames(pvals) = colnames(freq1)\n",
    "# write.csv(pvals,'healthy(control).csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in DR_data(counts):\n",
      "“not all rows sum up to 1 => normalization forced\n",
      "  some entries are 0 or 1 => transformation forced”\n"
     ]
    }
   ],
   "source": [
    "cov=read.csv('conv.csv',row.names=1)\n",
    "cov1 = data.frame(condition=factor(cov[rownames(freq),1], levels=c('mild(moderate)','healthy control','severe','convalescence','influenza')), \n",
    "                  row.names=rownames(freq))\n",
    "pvals = dirichlet_regression(counts=freq1, covariates=cov1, formula=counts ~ condition)$pvals\n",
    "colnames(pvals) = colnames(freq1)\n",
    "# write.csv(pvals,'mild(control).csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "heading_collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "R version 4.0.2 (2020-06-22)\n",
       "Platform: x86_64-pc-linux-gnu (64-bit)\n",
       "Running under: Ubuntu 16.04.6 LTS\n",
       "\n",
       "Matrix products: default\n",
       "BLAS:   /usr/lib/libblas/libblas.so.3.6.0\n",
       "LAPACK: /usr/lib/lapack/liblapack.so.3.6.0\n",
       "\n",
       "locale:\n",
       " [1] LC_CTYPE=C.UTF-8           LC_NUMERIC=C              \n",
       " [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    \n",
       " [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   \n",
       " [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 \n",
       " [9] LC_ADDRESS=C               LC_TELEPHONE=C            \n",
       "[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       \n",
       "\n",
       "attached base packages:\n",
       "[1] stats     graphics  grDevices utils     datasets  methods   base     \n",
       "\n",
       "other attached packages:\n",
       " [1] cowplot_1.1.1      forcats_0.5.1      stringr_1.4.0      dplyr_1.0.4       \n",
       " [5] purrr_0.3.4        readr_1.4.0        tidyr_1.1.2        tibble_3.0.6      \n",
       " [9] ggplot2_3.3.3      tidyverse_1.3.0    data.table_1.13.6  DirichletReg_0.7-0\n",
       "[13] Formula_1.2-4      Matrix_1.2-18      RColorBrewer_1.1-2 SeuratObject_4.0.0\n",
       "[17] Seurat_4.0.0      \n",
       "\n",
       "loaded via a namespace (and not attached):\n",
       "  [1] Rtsne_0.15           colorspace_2.0-0     deldir_0.2-10       \n",
       "  [4] ellipsis_0.3.1       ggridges_0.5.3       IRdisplay_1.0       \n",
       "  [7] fs_1.5.0             base64enc_0.1-3      rstudioapi_0.13     \n",
       " [10] spatstat.data_2.0-0  leiden_0.3.7         listenv_0.8.0       \n",
       " [13] ggrepel_0.9.1        lubridate_1.7.9.2    xml2_1.3.2          \n",
       " [16] codetools_0.2-16     splines_4.0.2        polyclip_1.10-0     \n",
       " [19] IRkernel_1.1.1       jsonlite_1.7.2       broom_0.7.5         \n",
       " [22] ica_1.0-2            dbplyr_2.1.0         cluster_2.1.0       \n",
       " [25] png_0.1-7            uwot_0.1.10          shiny_1.6.0         \n",
       " [28] sctransform_0.3.2    compiler_4.0.2       httr_1.4.2          \n",
       " [31] backports_1.2.1      assertthat_0.2.1     fastmap_1.1.0       \n",
       " [34] lazyeval_0.2.2       cli_2.3.0            later_1.1.0.1       \n",
       " [37] htmltools_0.5.1.1    tools_4.0.2          igraph_1.2.6        \n",
       " [40] gtable_0.3.0         glue_1.4.2           RANN_2.6.1          \n",
       " [43] reshape2_1.4.4       Rcpp_1.0.6           spatstat_1.64-1     \n",
       " [46] scattermore_0.7      cellranger_1.1.0     vctrs_0.3.6         \n",
       " [49] nlme_3.1-149         lmtest_0.9-38        ps_1.5.0            \n",
       " [52] globals_0.14.0       rvest_0.3.6          mime_0.10           \n",
       " [55] miniUI_0.1.1.1       lifecycle_1.0.0      irlba_2.3.3         \n",
       " [58] goftest_1.2-2        future_1.21.0        MASS_7.3-53         \n",
       " [61] zoo_1.8-8            scales_1.1.1         hms_1.0.0           \n",
       " [64] miscTools_0.6-26     promises_1.2.0.1     spatstat.utils_2.0-0\n",
       " [67] parallel_4.0.2       sandwich_3.0-0       reticulate_1.18     \n",
       " [70] pbapply_1.4-3        gridExtra_2.3        rpart_4.1-15        \n",
       " [73] stringi_1.5.3        repr_1.1.3           rlang_0.4.10        \n",
       " [76] pkgconfig_2.0.3      matrixStats_0.58.0   evaluate_0.14       \n",
       " [79] lattice_0.20-41      ROCR_1.0-11          tensor_1.5          \n",
       " [82] patchwork_1.1.1      htmlwidgets_1.5.3    tidyselect_1.1.0    \n",
       " [85] parallelly_1.23.0    RcppAnnoy_0.0.18     plyr_1.8.6          \n",
       " [88] magrittr_2.0.1       R6_2.5.0             generics_0.1.0      \n",
       " [91] pbdZMQ_0.3-5         DBI_1.1.1            withr_2.4.1         \n",
       " [94] haven_2.3.1          pillar_1.4.7         mgcv_1.8-33         \n",
       " [97] fitdistrplus_1.1-3   survival_3.2-3       abind_1.4-5         \n",
       "[100] future.apply_1.7.0   modelr_0.1.8         crayon_1.4.1        \n",
       "[103] uuid_0.1-4           KernSmooth_2.23-17   plotly_4.9.3        \n",
       "[106] maxLik_1.4-6         readxl_1.3.1         grid_4.0.2          \n",
       "[109] reprex_1.0.0         digest_0.6.27        xtable_1.8-4        \n",
       "[112] httpuv_1.5.5         munsell_0.5.0        viridisLite_0.3.0   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sessionInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
